function [s, Lambda] = omp(Phi, v, m)

% s = omp(Phi, v, m)
%
% Model:
% v = Phi * s
%
% Input argument:
% Phi is a N by d sensing matrix.
% v is an observed n-dimensional signal
% m is m-sparse
%
% Output argument:
% s is recoveried signal
% Lambda is a set of none zero location of s.
%  
% Reference:
% <<Signal Recovery From Random Measurements Via Orthogonal Matching Pursuit>>, Joel A. Tropp, and Anna c. Gilbert
%
% Author: 
% Cao Jin, xlscaoj@gmail.com
% 
% History:
% 2012/5/21: Initail release


% Initialize

  [N, d] = size(Phi);
  s_hat = zeros(d, 1);
  Lambda = [];
  r = v;
  Phit = [];

% Procedure

  for t = 1 : m
    % 
    [~, l] = max(abs(Phi' * r));

    % A set collects locations of signal s's none zero value
    Lambda = [Lambda l];

    % Phit should be the Moore-Penrose pseudo-inverse of Phi, the
    % following line should be correct!
    Phit = [Phit Phi(:, l)];

    [x, ~] = lsqlin(Phit, v);

    % 
    s_hat(l) = x(length(x));

    % update the resdual
    r = v - Phit * x ;
  end

s = s_hat;
